// test.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

#include <engine.h>
#include "math.h"
#include <stdlib.h>
#include <stdio.h>
#define PI 3.14159265
const int Idim = 1000;
const double teta = 0.05;

typedef double fieldComponent[Idim];
struct particle
    {
        double Xx;
        double Xy;
        double Xxb;
        double Xyb;
        double ux;
        double uy;
        double q;
        double m;
        int active;
        
    };
class fileout2 {
public:
  fileout2() { fp=fopen("sample.out","w+");}
  ~fileout2() { fclose(fp);}
  FILE *fp;
};

double min (double x, double y) {
    if (x>=y) return y; else return x;
}

double max (double x, double y) {
    if (x>=y) return x; else return y;
}

double fRand(double fMin, double fMax)
{
    double f = (double)rand() / RAND_MAX;
    return fMin + f * (fMax - fMin);
}

double trunc(double d){ return (d>0) ? floor(d) : ceil(d) ; }

void initializeParticles (particle* particles,double plasmaDensityPerM,double clusterFactorIn,int particlesPerCell, double ddx, int numOfParticles)  
{
    
    printf ("total particles: %d, each is a cluster of %e \n",numOfParticles,clusterFactorIn);
    int particleSerial=0;
    double clusterFactor;
    double T=1000;
    int i,j,k;
    for (i=50; i<Idim-50; i++) 
        for (j=Idim/4; j<Idim/4*3; j++) 
            for (k=0; k<particlesPerCell; k++) {
                // put ion cluster
                clusterFactor = clusterFactorIn;//if (i<200) clusterFactor = clusterFactorIn*(i-49)/150;
        
                particles[particleSerial].Xx=i*ddx+ddx/2;//(((double) (particleSerial%(particlesPerCell) +1))/(particlesPerCell+1))*ddx;
                particles[particleSerial].Xy=j*ddx+ddx/2;//(((double) (particleSerial%(particlesPerCell) +1))/(particlesPerCell+1))*ddx;
                particles[particleSerial].ux=fRand (-40e-5,40e-5);
                particles[particleSerial].uy=fRand (-40e-5,40e-5);
                particles[particleSerial].q=clusterFactor* 1.6e-19;
                particles[particleSerial].m=clusterFactor*14*1.66e-27;
                particles[particleSerial].active=1;
               

                // put electron cluster
                //i*ddx+ddx/2;//
                //j*ddx+ddx/2;
                particles[particleSerial+1].Xx=i*ddx+ddx/2;//(((double) (particleSerial%(particlesPerCell) +1))/(particlesPerCell+1))*ddx;
                particles[particleSerial+1].Xy=j*ddx+ddx/2;//(((double) (particleSerial%(particlesPerCell) +1))/(particlesPerCell+1))*ddx;
                particles[particleSerial+1].ux=fRand (-40,40);
                particles[particleSerial+1].uy=fRand (-40,40);
                particles[particleSerial+1].q = clusterFactor*-1.6e-19;
                particles[particleSerial+1].m = clusterFactor*9.1e-31 ;
                particles[particleSerial+1].active=1; 
               
                particleSerial=particleSerial+2;
                 
                 }  
    printf  ("done init particles \n");
}
 

// mover
// PIC particle mover
          void mover (particle* particles, int numOfParticles, fieldComponent* Ex,fieldComponent* Ey,fieldComponent* Hz,double ddx ,double dt, double speedOfLight )
             
          {
          double magneticPermeability0 = 4.0 * 3.14 * 1.0e-7;//permeability of free space    
          int i;    
          for (i=0; i<numOfParticles; i++){
                
                 if (particles[i].active==0) continue;
                 particles[i].Xxb = particles[i].Xx;
                 particles[i].Xyb = particles[i].Xy;
                                               
                 double ii = floor (particles[i].Xx/ddx);
                 double jj = floor (particles[i].Xy/ddx);
                 
                
                 
                 double wx = ii+1-particles[i].Xx/ddx;
                 double wxp =1-wx;
                 double wy = jj+1-particles[i].Xy/ddx;
                 double wyp = 1-wy;
                
                 
                 double ExEquivalent = Ex[(int) (particles[i].Xx/ddx)][(int) (particles[i].Xy/ddx)]*wx*wy;
                 ExEquivalent = ExEquivalent + Ex[((int) (particles[i].Xx/ddx))+1][(int) (particles[i].Xy/ddx)]*wxp*wy;
                 ExEquivalent = ExEquivalent + Ex[(int) (particles[i].Xx/ddx)][((int) (particles[i].Xy/ddx))+1]*wx*wyp;
                 ExEquivalent = ExEquivalent + Ex[((int) (particles[i].Xx/ddx))+1][((int) (particles[i].Xy/ddx))+1]*wxp*wyp;
              
               
                 double EyEquivalent = Ey[(int) (particles[i].Xx/ddx)][(int) (particles[i].Xy/ddx)]*wx*wy;
                 EyEquivalent = EyEquivalent + Ey[((int) (particles[i].Xx/ddx))+1][(int) (particles[i].Xy/ddx)]*wxp*wy;
                 EyEquivalent = EyEquivalent + Ey[(int) (particles[i].Xx/ddx)][((int) (particles[i].Xy/ddx))+1]*wx*wyp;
                 EyEquivalent = EyEquivalent + Ey[((int) (particles[i].Xx/ddx))+1][((int) (particles[i].Xy/ddx))+1]*wxp*wyp;
                 
                 
                 double HzEquivalent = Hz[(int) (particles[i].Xx/ddx)][(int) (particles[i].Xy/ddx)]*wx*wy;
                 HzEquivalent = HzEquivalent + Hz[((int) (particles[i].Xx/ddx))+1][(int) (particles[i].Xy/ddx)]*wxp*wy;
                 HzEquivalent = HzEquivalent + Hz[(int) (particles[i].Xx/ddx)][((int) (particles[i].Xy/ddx))+1]*wx*wyp;
                 HzEquivalent = HzEquivalent + Hz[((int) (particles[i].Xx/ddx))+1][((int) (particles[i].Xy/ddx))+1]*wxp*wyp;
                 
                         
                 HzEquivalent = HzEquivalent*magneticPermeability0;
                 double em = particles[i].q/particles[i].m;
                 double uxminus = particles[i].ux + em*ExEquivalent*dt/2;
                 double uyminus = particles[i].uy + em*EyEquivalent*dt/2;
                 double gammaX = sqrt(1+uxminus*uxminus/speedOfLight/speedOfLight);
                 double gammaY = sqrt(1+uyminus*uyminus/speedOfLight/speedOfLight);
                 double uxzero = uxminus + em*uyminus/gammaY* HzEquivalent*dt/2;
                 double uyzero = uyminus - em*uxminus/gammaX* HzEquivalent*dt/2;
                 gammaX = sqrt(1+uxzero*uxzero/speedOfLight/speedOfLight);
                 gammaY = sqrt(1+uyzero*uyzero/speedOfLight/speedOfLight);
                 
                 double uxplus = uxminus + 2*em/(1+(em*dt/2)*(em*dt/2)*HzEquivalent*HzEquivalent)*(uyzero/gammaY*HzEquivalent*dt/2);
                 double uyplus = uyminus + 2*em/(1+(em*dt/2)*(em*dt/2)*HzEquivalent*HzEquivalent)*-(uxzero/gammaX*HzEquivalent*dt/2);
                 particles[i].ux = uxplus + em*ExEquivalent*dt/2;
                 particles[i].uy = uyplus + em*EyEquivalent*dt/2;
                 gammaX = sqrt(1+particles[i].ux*particles[i].ux/speedOfLight/speedOfLight);
                 gammaY = sqrt(1+particles[i].uy*particles[i].uy/speedOfLight/speedOfLight);
                 particles[i].Xx= particles[i].Xx + dt*particles[i].ux/gammaX;
                 particles[i].Xy= particles[i].Xy + dt*particles[i].uy/gammaY;
                 
                 //if ((particles[i].Xx/ddx<=50) & (particles[i].ux<0)) particles[i].ux = -0.5*particles[i].ux;
                 //if ((particles[i].Xx/ddx>=Idim-50) & (particles[i].ux>0)) particles[i].ux = -0.5*particles[i].ux;
                 //if ((particles[i].Xy/ddx<=Idim/4) & (particles[i].uy<0)) particles[i].uy = -0.5*particles[i].uy;
                 //if ((particles[i].Xy/ddx>=Idim/4*3) & (particles[i].uy>0)) particles[i].uy = -0.5*particles[i].uy;
                 
                 if (((particles[i].Xx/ddx)>=Idim-1) | (particles[i].Xy/ddx>=Idim-1) | (particles[i].Xx/ddx<=1) | (particles[i].Xy/ddx<=1)) particles[i].active=0;  
                 
             }
        }
          
             void currentDiposition (particle* particles, fieldComponent* Jx, fieldComponent* Jy, int numOfParticles, double ddx, double dt)
          {
          int i,j;
          // pic current diposition
          
            for (i = 0; i < Idim; i++) 
                for (j = 0; j < Idim; j++) {
                    Jx[i][j]=0;
                    Jy[i][j]=0;
                }
           
            for (i=0; i<numOfParticles; i++){
                if (particles[i].active==0) continue;
                double x1 = particles[i].Xxb;
                double x2 = particles[i].Xx;
                double y1 = particles[i].Xyb;
                double y2 = particles[i].Xy;
                int i1 = floor (x1/ddx);
                int i2 = floor (x2/ddx);
                int j1 = floor (y1/ddx);
                int j2 = floor (y2/ddx);
                double xr = min (min (i1*ddx, i2*ddx) + ddx, max (max (i1*ddx,i2*ddx), (x1+x2)/2));
                double yr = min (min (j1*ddx, j2*ddx) + ddx, max (max (j1*ddx,j2*ddx), (y1+y2)/2));
                double q = particles[i].q;
                double fx1 = q/dt*(xr-x1);
                double fx2 = q/dt*(x2-xr);
                double fy1 = q/dt*(yr-y1);
                double fy2 = q/dt*(y2-yr);
                double wx1 = (x1+xr)/2/ddx-i1;
                double wx2 = (xr+x2)/2/ddx-i2;
                double wy1 = (y1+yr)/2/ddx-j1;
                double wy2 = (yr+y2)/2/ddx-j2;
                
                Jx[i1][j1] = Jx[i1][j1] + fx1*(1-wy1)/ddx/ddx;
                Jx[i1][j1+1] = Jx[i1][j1+1] + fx1*wy1/ddx/ddx;
                Jx[i2][j2] = Jx[i2][j2]+ fx2*(1-wy2)/ddx/ddx;
                Jx[i2][j2+1] = Jx[i2][j2+1]+ fx2*wy2/ddx/ddx;
                
                Jy[i1][j1] = Jy[i1][j1] + fy1*(1-wx1)/ddx/ddx;
                Jy[i1+1][j1] = Jy[i1+1][j1] + fy1*wx1/ddx/ddx;
                Jy[i2][j2] = Jy[i2][j2]+ fy2*(1-wx2)/ddx/ddx;
                Jy[i2+1][j2] = Jy[i2+1][j2]+ fy2*wx2/ddx/ddx;
               
             }
          
           
          }

          
         


void mexFunction( )
{
	// set MATLAB engine
	Engine *m_pEngine;
	m_pEngine = engOpen (NULL);

   // simulation size definition & consts
   const int iterations =500; // number of iterations.
   const int Jdim=Idim;
   const int CppToMatlabConst = 20; // the relation between c++ and MATLAB iterations.
   const int nPml = 8; // thickness of PML
   const double speedOfLight = 2.99792458e8;  
   const double ddx = 100e-9; // cell size in nm.
   const double dt = ddx/(2*speedOfLight); // time step.
   const double gradingOrder = 3; // grading order for PML
   const int sampleX = 975; // sample point for FFT 
   const int sampleY = 500;  
   double magneticPermeability0 = 4.0 * PI * 1.0e-7;//permeability of free space
   double electricalPermittivity0 = 1.0 / (speedOfLight * speedOfLight * magneticPermeability0);       //permittivity of free space
   double temperture[iterations]; 
   double sample[iterations];
   
   
   // experiment setup
   const double pulseEnergy =4e-3; //1e-7; // in J
   const double sourceWaveLength = 1e-6; // source wavelength in m
   const double freq =speedOfLight/sourceWaveLength; // frequancy in Hz
   const double pulseWidth = 33e-15; // pulse width 
   const double t0 =pulseWidth*3; // 3; //pulse middle point
   const double w0 = 4e-6;
   const double lf = 100e-6;
   const double A0 =  sqrt (2 * pulseEnergy/(pulseWidth*PI*w0*w0*speedOfLight*electricalPermittivity0));
   printf ("A0 %e \n",A0);
   // pulse 2
   const double pulseEnergy2 =15e-3; //1e-7; // in J
   const double sourceWaveLength2 = 1e-6; // source wavelength in m
   const double freq2 =speedOfLight/sourceWaveLength2; // frequancy in Hz
   const double pulseWidth2 =6.6*5e-15; // pulse width 
   const double t02 =pulseWidth*3; //pulse middle point
   const double w02 = 4e-6;
   const double lf2 = 300e-6;
   const double A02 =  sqrt (2 * pulseEnergy2/(pulseWidth2*PI*w02*w02*speedOfLight*electricalPermittivity0));
   
    // definition of variables and arrays
    double temporary;
    
    typedef double fieldComponent[Idim];
        fieldComponent* phi = new fieldComponent[Jdim];
        fieldComponent* Ex = new fieldComponent[Jdim];
        fieldComponent* Exb = new fieldComponent[Jdim];
        fieldComponent* Exbb = new fieldComponent[Jdim];
        fieldComponent* Exbbb = new fieldComponent[Jdim];
        fieldComponent* Ey = new fieldComponent[Jdim];
        fieldComponent* Eyb = new fieldComponent[Jdim];
        fieldComponent* Eybb = new fieldComponent[Jdim];
        fieldComponent* Eybbb = new fieldComponent[Jdim];
        fieldComponent* Jy = new fieldComponent[Jdim];
        fieldComponent* Jx = new fieldComponent[Jdim];
        fieldComponent* Hz = new fieldComponent[Jdim];
        fieldComponent* Dx = new fieldComponent[Jdim];
        fieldComponent* Dy = new fieldComponent[Jdim];
        fieldComponent* Bz = new fieldComponent[Jdim];
        fieldComponent* chargeDensity = new fieldComponent[Jdim];
    
    int i ,j,k, t, particleI;
    
    //PIC
    
    int notActive=0;
    int particlesPerCell=1;
    int numOfParticles = (Idim-100)*(Jdim/2)*particlesPerCell*2;
    
    particle* particles = new particle[numOfParticles];
    double plasmaDensityPerM = 3.5e25;
    double clusterFactor = (Idim/2*ddx)*(Jdim*ddx)*plasmaDensityPerM/numOfParticles/2;
    initializeParticles (particles, plasmaDensityPerM,clusterFactor,particlesPerCell,ddx,numOfParticles);
    // initialization
    
    for (i=0; i<Idim; i++)
        for (j=0; j<Jdim; j++)
        {
            Ex[i][j]=0;
            Ey[i][j]=0; 
            Jx[i][j]=0;
            Jy[i][j]=0;
            Hz[i][j]=0; 
            Dx[i][j]=0;
            Dy[i][j]=0;
            Bz[i][j]=0;
            phi[i][j]=0;
            
            
        }
   
     //openArray (phi);
     //electricFieldFromPhi (phi, Ex, Ey,ddx);
   
    // PML coefficients
    typedef double pmlCoefficentX[Idim];
        pmlCoefficentX* c2x = new pmlCoefficentX[Jdim];
        pmlCoefficentX* c3x = new pmlCoefficentX[Jdim];
    typedef double pmlCoefficentXHigh;
        pmlCoefficentXHigh* c4x = new pmlCoefficentXHigh[Idim];
        pmlCoefficentXHigh* c5x = new pmlCoefficentXHigh[Idim];
    typedef double pmlCoefficentY[Idim];
        pmlCoefficentY* c2y = new pmlCoefficentY[Jdim];
        pmlCoefficentY* c3y = new pmlCoefficentY[Jdim];
    typedef double pmlCoefficentYHigh;
        pmlCoefficentYHigh* c4y = new pmlCoefficentYHigh[Jdim];
        pmlCoefficentYHigh* c5y = new pmlCoefficentYHigh[Jdim];
    typedef double pmlCoefficentH;
        pmlCoefficentH* d1z = new pmlCoefficentH[Idim];
        pmlCoefficentH* d2z = new pmlCoefficentH[Idim];
        pmlCoefficentH* d3z = new pmlCoefficentH[Jdim];
        pmlCoefficentH* d4z = new pmlCoefficentH[Jdim];

    // initialize coefficents
    double gradientK = 1.0;
    double impedance0 = sqrt(magneticPermeability0 / electricalPermittivity0);
    double reflectionCoefficient0 = 1.0e-7;
    double boundaryWidth = (double  ) nPml * ddx;    // width of PML region (in mm)
    double electricalConductivityMaximum = -log(reflectionCoefficient0) * (gradingOrder + 1.0) / (2.0 * impedance0 * boundaryWidth);
    double boundaryFactor = electricalConductivityMaximum / ( ddx * (pow(boundaryWidth,gradingOrder)) * (gradingOrder + 1));
    double gradientConductivity; // not initialized!
    double x;
    double c1x = dt / (electricalPermittivity0 * ddx) ;
    double c1y = c1x;
    double temp = (electricalPermittivity0 ) / dt ;
    for (i=0; i<Idim; i++)
        for (j=0; j<Jdim; j++)
        {
            c2x[i][j]=1; 
            c3x[i][j]=1./temp;
            c2y[i][j]=1;
            c3y[i][j]=1./temp;
            c4x[i]=temp;
            c5x[i]=-temp;
            c4y[i]=temp;
            c5y[i]=-temp;
            d1z[i]=1;
            d2z[i]=dt/ddx;
            d3z[i]=1;
            d4z[i]=1.0 / magneticPermeability0;
        }
        
    // calculate gradient
    double gradientC2[nPml];
    double gradientC3[nPml];
    double gradientC4[nPml];
    double gradientC5[nPml];
    double gradientD1[nPml];
    double gradientD2[nPml];
    double gradientD3[nPml];
    double gradientD4[nPml];
  
    for (i = 0, x = 0.0; i < nPml; i++, x++) {
        double x1 = (x + 0.5) * ddx;       
        double x2 = (x - 0.5) * ddx;       
        if (i == 0) 
            gradientConductivity = boundaryFactor * (pow(x1,(gradingOrder+1))  );   //   polynomial grading  (special case: on the edge, 1/2 = pml, 1/2 = vacuum)
            else gradientConductivity = boundaryFactor * (pow(x1,(gradingOrder+1)) - pow(x2,(gradingOrder+1)) );   //   polynomial grading
        double temporary = gradientK * electricalPermittivity0 / dt ;
        double temporary1 = gradientConductivity / 2.0;
        gradientC2[i] = (temporary - temporary1) / (temporary + temporary1);
        gradientC3[i] = 1.0 / (temporary + temporary1);
        x1 = (x + 1.0) * ddx;       // upper bounds for point i
        x2 = (x + 0.0) * ddx;       // lower bounds for point i
        gradientConductivity = boundaryFactor * (pow(x1,(gradingOrder+1)) - pow(x2,(gradingOrder+1)) );   //   polynomial grading
        temporary = gradientK * electricalPermittivity0 / dt ;
        temporary1 = gradientConductivity / 2.0;
        gradientC4[i] = (temporary + temporary1);
        gradientC5[i] = (temporary1 - temporary);
        temporary = gradientK / dt ;
        temporary1 = gradientConductivity / (2.0 * electricalPermittivity0);
        gradientD1[i] = (temporary - temporary1) / (temporary + temporary1);
        gradientD2[i] = (1.0 / ddx) / (temporary + temporary1);
        temporary  *= magneticPermeability0;
        temporary1 *= magneticPermeability0;
        gradientD3[i] = (temporary - temporary1) / (temporary + temporary1);
        gradientD4[i] = (1.0 / dt) / (temporary + temporary1);
    } 

        int abcSize = nPml;  // for compatibility with PML code
        int xSize = Idim;
        int ySize = Jdim;
        for (j = 0; j < abcSize; j++) {                         
        

        for (i = 0; i < xSize; i++) {
            c2x[i][abcSize - j] = gradientC2[j];     // front
            c3x[i][abcSize - j] = gradientC3[j];
            c2x[i][ySize - abcSize + j] = gradientC2[j];     // back
            c3x[i][ySize - abcSize + j] = gradientC3[j];
        } 
        c4y[abcSize - j -1] = gradientC4[j];     // front
        c5y[abcSize - j -1] = gradientC5[j];
        c4y[ySize - abcSize + j] = gradientC4[j];     // back
        c5y[ySize - abcSize + j] = gradientC5[j];
        d3z[abcSize - j - 1] = gradientD3[j];     // front, 
        d4z[abcSize - j - 1] = gradientD4[j];     //   "         
        d3z[ySize - abcSize + j] = gradientD3[j];     // back
        d4z[ySize - abcSize + j] = gradientD4[j];     //   "
    } 
    // ey,hz --- left/right (x)
    for (i = 0; i < abcSize; i++) {                            // do coefficients for ey and hz
        // do coefficients for ey for left and right regions
        for (j = 0; j < ySize; j++) {
            c2y[abcSize - i][j] = gradientC2[i];     // left
            c3y[abcSize - i][j] = gradientC3[i];
            c2y[xSize - abcSize + i][j] = gradientC2[i];     // right
            c3y[xSize - abcSize + i][j] = gradientC3[i];
        }     
    
        c4x[abcSize - i -1] = gradientC4[i];     // left
        c5x[abcSize - i -1] = gradientC5[i];
        c4x[xSize - abcSize + i] = gradientC4[i];     // right
        c5x[xSize - abcSize + i] = gradientC5[i];
    
        d1z[abcSize - i - 1] = gradientD1[i];     // left, (note that the index is offset by 1 from c1y,c2y)
        d2z[abcSize - i - 1] = gradientD2[i];     //   "         
        d1z[xSize - abcSize + i] = gradientD1[i];     // right
        d2z[xSize - abcSize + i] = gradientD2[i];     //   "
    } 

    // temporary vars for pic
    double x1,x2,y1,y2,i1,i2,j1,j2,xr,yr,q,fx1,fx2,fy1,fy2,wx1,wx2,wy1,wy2;
    double gammaX1,gammaX2,gammaY1,gammaY2,tempux,hx,hy,ExEquivalent,EyEquivalent,HzEquivalent ;
    

    // main FDTD loop
    for (t=0; t<iterations; t++) 
    {
        
        // calculate electric fields
         // at upper pml
         printf ("iteration t %d \n",t);
         for (i = 1; i < Idim; i++) {
            for (j = 1; j < abcSize; j++) {       
                temporary = Dx[i][j];
                Dx[i][j] = Dx[i][j]+  c1x * ( Hz[i][j] - Hz[i][j-1] );
                Ex[i][j] = c2x[i][j] * Ex[i][j] + c3x[i][j] * (c4x[i] * Dx[i][j] + c5x[i] * temporary );
                temporary = Dy[i][j];
                Dy[i][j] = Dy[i][j]+ c1y * ( Hz[i-1][j] - Hz[i][j] );
                Ey[i][j] = c2y[i][j] * Ey[i][j] + c3y[i][j] * (c4y[j] * Dy[i][j] + c5y[j] * temporary );
            } 
        }
         // at lower pml
         for (i = 1; i < Idim; i++) {
            for (j = Jdim-abcSize; j < Jdim; j++) {       
                temporary = Dx[i][j];
                Dx[i][j] = Dx[i][j]+  c1x * ( Hz[i][j] - Hz[i][j-1] );
                Ex[i][j] = c2x[i][j] * Ex[i][j] + c3x[i][j] * (c4x[i] * Dx[i][j] + c5x[i] * temporary );
                temporary = Dy[i][j];
                Dy[i][j] = Dy[i][j]+ c1y * ( Hz[i-1][j] - Hz[i][j] );
                Ey[i][j] = c2y[i][j] * Ey[i][j] + c3y[i][j] * (c4y[j] * Dy[i][j] + c5y[j] * temporary );
            } 
        }
         // at left pml
         for (i = 1; i < abcSize; i++) {
            for (j = abcSize; j < Jdim-abcSize; j++) {       
                temporary = Dx[i][j];
                Dx[i][j] = Dx[i][j]+  c1x * ( Hz[i][j] - Hz[i][j-1] );
                Ex[i][j] = c2x[i][j] * Ex[i][j] + c3x[i][j] * (c4x[i] * Dx[i][j] + c5x[i] * temporary );
                temporary = Dy[i][j];
                Dy[i][j] = Dy[i][j]+ c1y * ( Hz[i-1][j] - Hz[i][j] );
                Ey[i][j] = c2y[i][j] * Ey[i][j] + c3y[i][j] * (c4y[j] * Dy[i][j] + c5y[j] * temporary );
            } 
        }
           // at right pml
         for (i = Idim-abcSize; i < Idim; i++) {
            for (j = abcSize; j < Jdim-abcSize; j++) {       
                temporary = Dx[i][j];
                Dx[i][j] = Dx[i][j]+  c1x * ( Hz[i][j] - Hz[i][j-1] );
                Ex[i][j] = c2x[i][j] * Ex[i][j] + c3x[i][j] * (c4x[i] * Dx[i][j] + c5x[i] * temporary );
                temporary = Dy[i][j];
                Dy[i][j] = Dy[i][j]+ c1y * ( Hz[i-1][j] - Hz[i][j] );
                Ey[i][j] = c2y[i][j] * Ey[i][j] + c3y[i][j] * (c4y[j] * Dy[i][j] + c5y[j] * temporary );
            } 
        }
            
            // center
         for (i = abcSize; i < Idim-abcSize; i++) {
            for (j = abcSize; j < Jdim-abcSize; j++) {     
                Exbbb[i][j] =Exbb[i][j];
                Exbb[i][j] =Exb[i][j];
                Exb[i][j] = Ex[i][j];
                Eybbb[i][j] =Eybb[i][j];
                Eybb[i][j] =Eyb[i][j];
                Eyb[i][j] = Ey[i][j];
                Ex[i][j] = Ex[i][j] +c1x* ( Hz[i][j] - Hz[i][j-1] )- Jx[i][j]*dt/electricalPermittivity0; 
                Ey[i][j] = Ey[i][j] +c1y* ( Hz[i-1][j] - Hz[i][j])- Jy[i][j]*dt/electricalPermittivity0;
            } 
        } 
           
        // correction due to TFSF to the electric field;
         
         for (j=0; j<Jdim; j++ ) 
         {
             double x = ddx * j;
             double x0 = Jdim/2 * ddx;
             double td = -lf / speedOfLight * (sqrt (1.+(x-x0)*(x-x0)/(lf*lf))-1.);
             Ey[nPml+1][j] = Ey[nPml+1][j] +  A0*exp (- (t*dt- t0) * (t*dt- t0)/ (2*pulseWidth*pulseWidth) )*sin (2*PI*freq*(dt*t-td))* exp (-(x-x0)*(x-x0)/(2*w0*w0));
             // second pulse
             //double td2 = -lf2 / speedOfLight * (sqrt (1.+(x-x0)*(x-x0)/(lf2*lf2))-1.);
             //Ey[Jdim - nPml-1][j] = Ey[Jdim - nPml-1][j] +  A02*exp (- (t*dt- t02) * (t*dt- t02)/ (2*pulseWidth2*pulseWidth2) )*sin (2*PI*freq2*(dt*t-td2))* exp (-(x-x0)*(x-x0)/(2*w02*w02));
         }
      
        // calculate magnetic field
       // at upper PML  
           for (i = 0; i < Idim; i++) {
            for (j = 0; j < abcSize-1; j++) {       
              temporary = Bz[i][j];
                Bz[i][j] = d1z[i] * Bz[i][j] + d2z[i] * ( Ex[i][j+1] - Ex[i][j] + Ey[i][j] - Ey[i+1][j] );
                Hz[i][j] = d3z[j] * Hz[i][j] + d4z[j] * ( Bz[i][j] - temporary );
            } 
        }
         // at lower pml
         for (i = 0; i < Idim-1; i++) {
            for (j = Jdim-abcSize; j < Jdim-1; j++) {       
                temporary = Bz[i][j];
                Bz[i][j] = d1z[i] * Bz[i][j] + d2z[i] * ( Ex[i][j+1] - Ex[i][j] + Ey[i][j] - Ey[i+1][j] );
                Hz[i][j] = d3z[j] * Hz[i][j] + d4z[j] * ( Bz[i][j] - temporary );
            } 
        }
         // at left pml
         for (i = 0; i < abcSize; i++) {
            for (j = abcSize; j < Jdim-abcSize; j++) {       
               temporary = Bz[i][j];
                Bz[i][j] = d1z[i] * Bz[i][j] + d2z[i] * ( Ex[i][j+1] - Ex[i][j] + Ey[i][j] - Ey[i+1][j] );
                Hz[i][j] = d3z[j] * Hz[i][j] + d4z[j] * ( Bz[i][j] - temporary );
            } 
        }
           // at right pml
         for (i = Idim-abcSize; i < Idim-1; i++) {
            for (j = abcSize; j < Jdim-abcSize; j++) {       
               temporary = Bz[i][j];
                Bz[i][j] = d1z[i] * Bz[i][j] + d2z[i] * ( Ex[i][j+1] - Ex[i][j] + Ey[i][j] - Ey[i+1][j] );
                Hz[i][j] = d3z[j] * Hz[i][j] + d4z[j] * ( Bz[i][j] - temporary );
            } 
        }
            
            // center
      for (i = abcSize; i < Idim-abcSize; i++) {
            for (j = abcSize; j < Jdim-abcSize; j++) {     
               temporary = Bz[i][j];
               double temp = (1+teta/2)*(Ex[i][j+1] - Ex[i][j] + Ey[i][j] - Ey[i+1][j]) - teta*(1-teta/2)*(Exb[i][j+1] - Exb[i][j] + Eyb[i][j] - Eyb[i+1][j]);
               temp = temp + 0.5*(1-teta)*(1-teta)*teta*(Exbb[i][j+1] - Exbb[i][j] + Eybb[i][j] - Eybb[i+1][j]);
               temp = temp + 0.5*(1-teta)*(1-teta)*teta*teta*(Exbbb[i][j+1] - Exbbb[i][j] + Eybbb[i][j] - Eybbb[i+1][j]);
               
                Hz[i][j] = Hz[i][j] + ( temp )*dt/ddx/magneticPermeability0;
              
            } 
        } 
         
       
         
         // move particles
         
         //mover (particles, numOfParticles, Ex, Ey, Hz,ddx,dt, speedOfLight);
         //currentDiposition (particles, Jx, Jy, numOfParticles, ddx, dt);    
     
        // sample electric field
         sample[t] = Ey [sampleX][sampleY];
         
        // send data to MATLAB for plotting
               
        if (t%CppToMatlabConst == 0) {
             
              /*  
             
            // plot 2d
                
            mxArray* ExToMat;
            int size[2] = {Idim, Jdim};
            ExToMat = mxCreateNumericArray(2,(mwSize*) size, mxDOUBLE_CLASS, mxREAL );
            double *tempd = mxGetPr(ExToMat);
            int posCounter=0;
            for (i=0; i<Idim; i++)
                for (j=0; j<Jdim; j++) 
                {
                     tempd[posCounter]= Ey[i][j];
                     posCounter++;
                }
            mxArray *lhs[1];
            lhs [0] = ExToMat;
            engPutVariable (eng, "x", ExToMat);
			engEvalString (eng,"plot2d (x);");
			//mexCallMATLAB(0, NULL, 1, lhs, "plot2d");
            mxDestroyArray  (ExToMat);
            printf ("Iteration number: %d \n",t);
            
            
            // plot electron density
            
            for (i=0; i<Idim; i++) for (j=0; j<Jdim; j++) chargeDensity[i][j]=0;
            
            for (i=0; i<numOfParticles; i++) {
                if (particles[i].q>0) continue;
                if (particles[i].active==0) continue;
                double hx = particles[i].Xx/ddx - floor(particles[i].Xx/ddx);
                double hy = particles[i].Xy/ddx - floor(particles[i].Xy/ddx);
                chargeDensity[(int) (particles[i].Xx/ddx)][(int) (particles[i].Xy/ddx)]=chargeDensity[(int) (particles[i].Xx/ddx)][(int) (particles[i].Xy/ddx)]+(1-hx)*(1-hy)*particles[i].ux/speedOfLight;
                chargeDensity[((int) (particles[i].Xx/ddx))+1][(int) (particles[i].Xy/ddx)]=chargeDensity[((int) (particles[i].Xx/ddx))+1][(int) (particles[i].Xy/ddx)]+hx*(1-hy)*particles[i].ux/speedOfLight;
                chargeDensity[(int) (particles[i].Xx/ddx)][((int) (particles[i].Xy/ddx))+1]=chargeDensity[(int) (particles[i].Xx/ddx)][((int) (particles[i].Xy/ddx))+1]+(1-hx)*hy*particles[i].ux/speedOfLight;
                chargeDensity[((int) (particles[i].Xx/ddx))+1][((int) (particles[i].Xy/ddx))+1]=chargeDensity[((int) (particles[i].Xx/ddx))+1][((int) (particles[i].Xy/ddx))+1]+hx*hy*particles[i].ux/speedOfLight;
                 
                
            }
           
            mxArray* ExToMat3;
            size[0] = Idim;
            size[1] = Jdim;
            ExToMat3 = mxCreateNumericArray(2, (mwSize*) size, mxDOUBLE_CLASS, mxREAL );
            tempd = mxGetPr(ExToMat3);
            posCounter=0;
            for (i=0; i<Idim; i++)
                for (j=0; j<Jdim; j++) 
                {
                     
                     tempd[posCounter]= chargeDensity[i][j];
                     posCounter++;
                }
            engPutVariable (eng, "x", ExToMat3);
			engEvalString (eng,"plotParticle (x);");
            lhs [0] = ExToMat3;
            //mexCallMATLAB(0, NULL, 1, lhs, "plotParticle");
            mxDestroyArray  (ExToMat3);
            */
        
        }
        
         
      
        
       
    
  
   

        
        
           
        
    } // end of main loop

    // close field sample file
   
  
   
    fileout2 fieldSample;
    for (i=0; i<iterations; i++ ) fprintf(fieldSample.fp,"%e \n", sample[i]);
    fclose(fieldSample.fp);
    
   
    // clean up heap
    
    delete [] Ex;
    delete [] Exb;
    delete [] Exbb;
    delete [] Exbbb;
    delete [] Ey;
    delete [] Eyb;
    delete [] Eybb;
    delete [] Eybbb;
    delete [] Jx;
    delete [] Jy;
  
    delete [] Hz;
    delete [] Dx;
    delete [] Dy;
    delete [] Bz;
    delete [] c2x;
    delete [] c3x;
    delete [] c4x;
    delete [] c5x;
    delete [] c2y;
    delete [] c3y;
    delete [] c4y;
    delete [] c5y;
    delete [] d1z;
    delete [] d2z;
    delete [] d3z;
    delete [] d4z;
    delete [] phi;
    delete [] chargeDensity;
    
    return;
 }

int _tmain(int argc, _TCHAR* argv[])
{
	printf ("bye world\n");
	mexFunction();
	printf ("bye world\n");
	return 0;
}

